* load the replication data set BCG_JOP_ReplicationData.dta, but do not save over it	

cap which clarify.hlp
if _rc ssc install clarify

	
capture program drop nnsim
program define nnsim, rclass
        version 16.1
        syntax [, yy(integer 1946)]
		tempvar temp0 temp1 temp2 temp3
        gen `temp0' = rbinomial(1, proatt0) if year_a ==`yy'
        summarize `temp0'
        return scalar att0 = r(sum)
        gen `temp1' = rbinomial(1, proatt1) if year_a ==`yy'
        summarize `temp1'
        return scalar att1 = r(sum)
        gen `temp2' = rbinomial(1, protol0) if year_a ==`yy'
        summarize `temp2'
        return scalar tol0 = r(sum)
        gen `temp3' = rbinomial(1, protol1) if year_a ==`yy'
        summarize `temp3'
        return scalar tol1 = r(sum)
    end



xtset audience year, yearly

	
	
capture drop b1-b4

logit audience_accel  Nat_total_full_5 tol_end_total_full_5 dll_total_full_5  if aec2==1 & year_a>1945 & (audience_pursue ==0 | pursue_start==1), vce(cluster audience)


keep if e(sample)

capture drop proatt*
qui gen proatt0=.
qui gen proatt1=.

estsimp logit audience_accel  Nat_total_full_5 tol_end_total_full_5 dll_total_full_5  if aec2==1 & year_a>1945 & (audience_pursue ==0 | pursue_start==1), vce(cluster audience)

forval oo = 1/5681 {
capture drop temp0 temp1
setx [`oo'] 
*setx at_total_full_5 0 tol_end_total_full_5 0 dll_total_full_5 0 
qui simqi, genpr(temp0 temp1)
qui sum temp1
qui replace proatt0 = r(mean) in `oo'


capture drop temp0 temp1
setx [`oo'] 
setx Nat_total_full_5 Nat_total_full_5[`oo']+1 
qui simqi, genpr(temp0 temp1)
qui sum temp1
qui replace proatt1 = r(mean) in `oo'

}



capture drop b1-b4

capture drop protol*

qui gen protol0=.
qui gen protol1=.

estsimp logit audience_accel  Nat_total_full_5 tol_end_total_full_5 dll_total_full_5  if aec2==1 & year_a>1945 & (audience_pursue ==0 | pursue_start==1), vce(cluster audience)

forval oo = 1/5681 {
capture drop temp0 temp1
setx [`oo'] 
qui simqi, genpr(temp0 temp1)
qui sum temp1
qui replace protol0 = r(mean) in `oo'


capture drop temp0 temp1
setx [`oo'] 
setx tol_end_total_full_5 tol_end_total_full_5[`oo']+1 
qui simqi, genpr(temp0 temp1)
qui sum temp1
qui replace protol1 = r(mean) in `oo'

}


qui gen proattdif = proatt1-proatt0
qui gen proattperc = (proatt1-proatt0)/proatt0
qui gen protoldif = protol1-protol0
qui gen protolperc = (protol1-protol0)/protol0



capture drop Natt* Ntol*

qui gen Natt0=.
qui gen Natt1=.
qui gen Ntol0=.
qui gen Ntol1=.

qui gen Natt0m=.
qui gen Natt1m=.
qui gen Ntol0m=.
qui gen Ntol1m=.


forval yyy = 1946/2018 {
di `yyy'
capture drop simtemtol0 simtemtol1 simtematt0 simtematt1 
qui gen simtemtol0=.
qui gen simtemtol1=.
qui gen simtematt0=.
qui gen simtematt1=.
forval it = 1/1000 {
qui nnsim, yy(`yyy')	
qui replace simtemtol0 = r(tol0) in `it'
qui replace simtemtol1 = r(tol1) in `it'
qui replace simtematt0 = r(att0) in `it'
qui replace simtematt1 = r(att1) in `it'
}
qui sum simtemtol0
qui replace Ntol0= r(mean) if year_a == `yyy'
qui sum simtemtol1
qui replace Ntol1= r(mean) if year_a == `yyy'
qui sum simtematt0
qui replace Natt0= r(mean) if year_a == `yyy'
qui sum simtematt1
qui replace Natt1= r(mean) if year_a == `yyy'

qui centile simtemtol0
qui replace Ntol0m= r(c_1) if year_a == `yyy'
qui centile simtemtol1
qui replace Ntol1m= r(c_1) if year_a == `yyy'
qui centile simtematt0
qui replace Natt0m= r(c_1) if year_a == `yyy'
qui centile simtematt1
qui replace Natt1m= r(c_1) if year_a == `yyy'

}



keep audience year_a audience_on audience_explore audience_pursue audience_start audience_end  Nat_total_full_5 tol_end_total_full_5 dll_total_full_5 aec2 proatt0-Ntol1m


duplicates drop year_a, force

tsset year_a
tssmooth ma Natt05= Natt0 , window(5 1)
replace Natt05 = Natt05*6

tssmooth ma Natt15= Natt1 , window(5 1)
replace Natt15 = Natt15*6


tssmooth ma Ntol05= Ntol0 , window(5 1)
replace Ntol05 = Ntol05*6

tssmooth ma Ntol15= Ntol1 , window(5 1)
replace Ntol15 = Ntol15*6

*save "Figure1Output.dta", replace

twoway (line Ntol05 year_a)(line Ntol15 year_a), legend(label(1 "Observed Expectation") label(2 "With an Extra Toleration")) xtitle("Year") ytitle("Expected no of States with Increased Activity")

*graph save "Graph" "Figure1a.gph"

twoway (line Natt05 year_a)(line Natt15 year_a), legend(label(1 "Observed Expectation") label(2 "With an Extra Attack")) xtitle("Year") ytitle("Expected no of States with Increased Activity")

*graph save "Graph" "Figure1b.gph"
